#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <GL/glew.h>
#include <GL/glut.h>
#include <GL/gl.h>
#include <GL/glu.h>
#include <math.h>
#define pi 3.1415926535897
#include <GL/glxew.h>
#include <GL/glx.h>
#include <unistd.h>
#include <sys/time.h>

// global vars //
static int winWidth, winHeight;
static double winRatio;


const GLfloat black[] = { 0.0f, 0.0f, 0.0f, 1.0f };
const GLfloat white[] = { 1.0f, 1.0f, 1.0f, 1.0f };

static const int dataSize = 512;
static float* data;
static GLuint transform[3], spectrum, palette;
static int currentTransform = 0;
static bool ready = true;
static bool fast = true;
static float fading = 0.0;
static GLuint mapper;

// Utilities && Low level //

inline float abs(float x) {return x < 0 ? -x : x; }
inline double abs(double x) {return x < 0 ? -x : x; }

// standard texture setup
void textureSetup() {
	glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
	glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
	glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);	
}

// round to nearest power of 2
int round2(int x) {
	int n = 0;
	while(x >>= 1) n++;
	return 1 << n;
}

// time counter
class Chrono {
    private:
        double start_time;
    public:
        static double get_time() {
            struct timeval tp;
            struct timezone tzp;
            gettimeofday(&tp, &tzp);
            return ((double)tp.tv_sec+(1.e-6)*tp.tv_usec);
        }
        void start(){start_time = get_time();}
        double time(){return get_time() - start_time;}
        Chrono(){start();}
} chrono;

// standard projection setups
void perspectiveProjection() {
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glFrustum(-winRatio, winRatio, -1.0, 1.0, 2.0, 100.0);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	glEnable(GL_DEPTH_TEST);
	glDepthFunc(GL_LEQUAL);
}

void pixelProjection() {
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(0, winWidth, 0, winHeight, 0.0, 10.0);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	glDisable(GL_DEPTH_TEST);
}

void noProjection() {
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(-1.0, 1.0, -1.0, 1.0, 0.0, 10.0);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	glDisable(GL_DEPTH_TEST);
}

// GPU program handling
bool compileShader(int shader) {
	int l;
	glCompileShader(shader);
	glGetShaderiv(shader, GL_COMPILE_STATUS, &l);
	if(!l) {
		char * log = new char[1024];
		glGetShaderInfoLog(shader, 1024, &l, log); log[l] = 0;
		printf("Shader did not compile -- log follows: \n%s\n", log);
		delete[] log;
		return false;
	}
	return true;
}

int createProgram(const char* vertexCode, const char* pixelCode) {
	int vs = 0, ps = 0, l;
	if(vertexCode) {
		vs = glCreateShader(GL_VERTEX_SHADER);
		glShaderSource(vs, 1, &vertexCode, 0);
		if(!compileShader(vs)) {
			glDeleteShader(vs);
			vs = 0;
		}
	}
	if(pixelCode) {
		ps = glCreateShader(GL_FRAGMENT_SHADER);
		glShaderSource(ps, 1, &pixelCode, 0);
		if(!compileShader(ps)) {
			glDeleteShader(ps);
			ps = 0;
		}
	}
	int program = 0;
	if(vs || ps) {
		program = glCreateProgram();
		if(vs) glAttachShader(program, vs);
		if(ps) glAttachShader(program, ps);
		glLinkProgram(program);
		glGetProgramiv(program, GL_LINK_STATUS, &l);
		if(!l) {
			char * log = new char[1024];
			glGetProgramInfoLog(program, 1024, &l, log); log[l] = 0;
			printf("Program did not link -- log follows: \n%s\n", log);
			delete[] log;
			glDeleteProgram(program);
			program = 0;
		}
	}
	return program;
}

bool installProgram(int prog) {
	int l;
	if(prog) {
		glValidateProgram(prog);
		glGetProgramiv(prog, GL_VALIDATE_STATUS, &l);
	} else l = 1;
	if(l) glUseProgram(prog);
	else {
		char * log = new char[1024];
		glGetProgramInfoLog(prog, 1024, &l, log); log[l] = 0;
		printf("Program %d did not validate -- log follows: \n%s\n", prog, log);
		delete[] log;
	}
	return (bool)l;
}

bool checkFbo() {
	GLenum status = glCheckFramebufferStatusEXT(GL_FRAMEBUFFER_EXT);
	if(status != GL_FRAMEBUFFER_COMPLETE_EXT) {
		switch(status) {
			case GL_FRAMEBUFFER_INCOMPLETE_ATTACHMENT_EXT: puts("fbo incomplete"); break;
			case GL_FRAMEBUFFER_INCOMPLETE_DIMENSIONS_EXT: puts("fbo size mismatch"); break;
			case GL_FRAMEBUFFER_INCOMPLETE_MISSING_ATTACHMENT_EXT: puts("fbo empty"); break;
			case GL_FRAMEBUFFER_UNSUPPORTED_EXT: puts("fbo unsuported settings"); break;
		}
		return false;
	} else return true;
}

// basic geometry presets //
void screen() {
	glBegin(GL_QUADS);
	glTexCoord2f(0, 1); glVertex2f(-1, 1);
	glTexCoord2f(0, 0); glVertex2f(-1, -1);
	glTexCoord2f(1, 0); glVertex2f(1, -1);
	glTexCoord2f(1, 1); glVertex2f(1, 1);
	glEnd();
}


// GPU optimized FFT of 2D square float maps //
class FFT {
	private:
		bool is2D;
		int size, logSize;
		GLuint trigmap, permap, tex[2];
		static GLuint scrambler, butterfly[2], lastStage[2], inverter, fbo;
		int currentStep;
		int currentDir;
		int currentTex;
		float currentLog;
		GLuint tex1(float* data, GLenum format, GLenum components) {
			GLuint t;
			glGenTextures(1, &t);
			glBindTexture(GL_TEXTURE_1D, t);
			glTexImage1D(GL_TEXTURE_1D, 0, format, size, 0, components, GL_FLOAT, data);
			glTexParameterf(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_REPEAT);
			glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
			glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
			glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);			
			return t;
		}
		GLuint tex2(float* data, GLenum format, GLenum components) {
			GLuint t;
			glGenTextures(1, &t);
			glBindTexture(GL_TEXTURE_2D, t);
			glTexImage2D(GL_TEXTURE_2D, 0, format, size, is2D ? size : 1, 0, components, GL_FLOAT, data);
			glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
			glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
			glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);			
			return t;
		}
		
		// init offscreen rendering
		void initialize() {
			glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, fbo);
			glViewport(0, 0, size, is2D ? size : 1);
			noProjection();
		}
		
		// restore normal rendering state
		void finalize() {
			glActiveTextureARB(GL_TEXTURE0_ARB);
			installProgram(0);
			glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);		
		}
		
		// initial scrambling
		void _begin(GLuint in, bool inv = false) {
			glBindTexture(GL_TEXTURE_2D, in);
			if(inv) {
				glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_2D, tex[1], 0);
				installProgram(inverter);
				glUniform1f(glGetUniformLocation(inverter, "inv"), 1.0 / (float)(is2D ? size * size : size));
				screen();
				glBindTexture(GL_TEXTURE_2D, tex[1]);
			}
			glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_2D, tex[0], 0);
			glActiveTextureARB(GL_TEXTURE1_ARB);
			glBindTexture(GL_TEXTURE_1D, permap);
			installProgram(scrambler);
			glUniform1i(glGetUniformLocation(scrambler, "input"), 0);
			glUniform1i(glGetUniformLocation(scrambler, "perm"), 1);
			screen();
			currentStep = 0;
			currentDir = 0;
			currentTex = 0;
		}
		
		// iterative butterfly stages
		bool _end(GLuint out) {
			bool r = false;
			if(!currentStep++) {
				currentLog = 1;
			}
			if(currentStep <= logSize - 1) {
				GLuint b = butterfly[currentDir];
				installProgram(b);
				glUniform1i(glGetUniformLocation(b, "input"), 0);
				glUniform1i(glGetUniformLocation(b, "trig"), 1);
				glUniform1f(glGetUniformLocation(b, "margin"), 0.5 / (float)size);
				glActiveTextureARB(GL_TEXTURE1_ARB);
				glBindTexture(GL_TEXTURE_1D, trigmap);				
				glUniform1f(glGetUniformLocation(b, "inv"), 1.0 / (currentLog * size));
				currentLog *= 0.5;
				glUniform1f(glGetUniformLocation(b, "scale"), currentLog * size);
				glActiveTextureARB(GL_TEXTURE0_ARB);
				glBindTexture(GL_TEXTURE_2D, tex[currentTex++ & 1]);
				glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_2D, tex[currentTex & 1], 0);
				screen();
			} else {
				if(is2D) {
					if(currentDir) r = true;
					else currentStep = 0;
				} else r = true;
				installProgram(lastStage[currentDir]);
				glBindTexture(GL_TEXTURE_2D, tex[currentTex++ & 1]);
				glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_2D, r ? out : tex[currentTex & 1], 0);
				screen();
				currentDir++;
			}
			return r;
		}

	public:
		FFT(int s, bool _is2D = false) : is2D(_is2D), size(s) {
			float* buf = new float[size * 4];
			float invSize = 1 / (float)size;
			for(logSize = -1; s; s >>= 1) logSize++;
			
			// trig map
			for(int i = 0; i < size; i++) {
				buf[i * 4    ] = cos(2 * i * pi * invSize);
				buf[i * 4 + 1] = sin(2 * i * pi * invSize);
			}
			trigmap = tex1(buf, GL_RGBA16F_ARB, GL_RGBA);
			
			// permutation map
			bool done[size];
			for(int i = 0; i < size; i++) done[i] = false;
			for(int i = 0; i < size; i++) if(!done[i]) {
				int j = 0;
				int t = i;
				int base = size >> 1;
				for(int p = 0; p < logSize; p++) {
					if(t & 1) j += base;
					base >>= 1;
					t >>= 1;
				}
				buf[i] = (j + 0.5) * invSize;
				buf[j] = (i + 0.5) * invSize;
				done[i] = done[j] = true;
			}
			permap = tex1(buf, GL_ALPHA16, GL_ALPHA);
			delete[] buf;

			// temp textures allocation
			for(int i = 0; i < 2; i++) tex[i] = tex2(0, GL_RGBA16F_ARB, GL_RGBA);
			
			// compile shaders
			if(!fbo) { // not thread safe !
				glGenFramebuffersEXT(1, &fbo);
				
				inverter = createProgram(0,
					"uniform sampler2D input;"
					"uniform float inv;"
					"void main() {"
					"gl_FragColor = vec4(texture2D(input, gl_TexCoord[0].st).rg * vec2(1, -1) * inv, 0, 1);"
					"}"
				);
				
				scrambler = createProgram(0,
					"uniform sampler2D input;"
					"uniform sampler1D perm;"
					"void main() {"
					"gl_FragColor = texture2D(input, vec2(texture1D(perm, gl_TexCoord[0].s).a, texture1D(perm, gl_TexCoord[0].t).a));"
					"}"
				);

				butterfly[0] = createProgram(0,
					"uniform sampler2D input;"
					"uniform sampler1D trig;"
					"uniform float scale, inv, margin;"
					"vec2 cm(vec2 a, vec2 b) { return vec2(a.r * b.r - a.g * b.g, a.r * b.g + a.g * b.r); }"
					"void main() {"
					"float src = gl_TexCoord[0].s, lin = gl_TexCoord[0].t;"
					"float sig = sign(0.5 - frac(src * scale));"
					"vec2 butterfly = texture2D(input, vec2(src, lin)) * sig + texture2D(input, vec2(src + sig * inv, lin)).rg;"
					"vec2 result = frac(0.5 * src * scale) < 0.5 ? butterfly : cm(butterfly, texture1D(trig, 0.5 * frac(scale * (src - margin))));"
					"gl_FragColor = vec4(result, 0, 1);"
					"}"
				);
				
				butterfly[1] = createProgram(0,
					"uniform sampler2D input;"
					"uniform sampler1D trig;"
					"uniform float scale, inv, margin;"
					"vec2 cm(vec2 a, vec2 b) { return vec2(a.r * b.r - a.g * b.g, a.r * b.g + a.g * b.r); }"
					"void main() {"
					"float src = gl_TexCoord[0].t, lin = gl_TexCoord[0].s;"
					"float sig = sign(0.5 - frac(src * scale));"
					"vec2 butterfly = texture2D(input, vec2(lin, src)) * sig + texture2D(input, vec2(lin, src + sig * inv)).rg;"
					"vec2 result = frac(0.5 * src * scale) < 0.5 ? butterfly : cm(butterfly, texture1D(trig, 0.5 * frac(scale * (src - margin))));"
					"gl_FragColor = vec4(result, 0, 1);"
					"}"
				);
				
				lastStage[0] = createProgram(0,
					"uniform sampler2D input;"
					"void main() {"
					"vec2 src = gl_TexCoord[0].st;"
					"gl_FragColor = src.s < 0.5 ? texture2D(input, src) + texture2D(input, src + vec2(0.5, 0)) : texture2D(input, src + vec2(-0.5, 0)) - texture2D(input, src);"
					"}"
				);
				lastStage[1] = createProgram(0,
					"uniform sampler2D input;"
					"void main() {"
					"vec2 src = gl_TexCoord[0].st;"
					"gl_FragColor = src.t < 0.5 ? texture2D(input, src) + texture2D(input, src + vec2(0, 0.5)) : texture2D(input, src + vec2(0, -0.5)) - texture2D(input, src);"
					"}"
				);

			}
		}
		
		void begin(GLuint in, bool inv = false) {
			initialize();
			_begin(in, inv);
			finalize();
		}
		bool end(GLuint out) {
			initialize();
			bool r = _end(out);
			finalize();
			return r;
		}
		
		// expect 2 RGBA (size x 1) or (size x size) textures
		void compute(GLuint in, GLuint out, bool inv = false) {
			initialize();
			for(_begin(in, inv);!_end(out););
			finalize();
		}
};
GLuint FFT::scrambler, FFT::butterfly[2], FFT::lastStage[2], FFT::inverter, FFT::fbo = 0;
FFT* fft;


// input events handling //
void key(unsigned char k, int x, int y) {
	switch(k) {
		case 27: case 'q': exit(0);
	}
}

void displayFunc() { static int F = 0; chrono.start();
	if(!ready && fft->end(transform[currentTransform])) ready = true;
	fading += fast ? 0.1 : 1e-3;
	if(fading > 1) {
		// finish transform
		while(!ready) if(fft->end(transform[currentTransform])) ready = true;
		ready = false;
		
		// update state
		fading -= 1;
		if(currentTransform == 2) {
			currentTransform = 0; fast = false;
		} else currentTransform++;
		
		double order = rand() /(float)RAND_MAX - 1.5;
		double norm = 2e-5 * (order + 1) / (pow(0.5, order + 1) - pow(1e-6, order + 1)) / -order;
		for(int y = 0; y < dataSize; y++) {
			double fy = 0.5 - abs(0.5 - y / (double)dataSize);
			for(int x = 0; x < dataSize; x++) {
				double fx = 0.5 - abs(x / (double)dataSize - 0.5);
				double r = fx * fx + fy * fy;
				double p = r - 0.1; p *= p;
				double I = r > 1e-6 ? norm * pow(r, order): 0;
				double theta = rand() /(float)RAND_MAX * 2*pi;
				data[(dataSize * y + x) * 4] = I * cos(theta);
				data[(dataSize * y + x) * 4 + 1] = I * sin(theta);
			}
		}
		glBindTexture(GL_TEXTURE_2D, spectrum);
		glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, dataSize, dataSize, GL_RGBA, GL_FLOAT, data);
		fft->begin(spectrum);
	}
	
	// screen rendering setup
	glClearColor(0, 0, 1, 1);
	glViewport(0, 0, winWidth, winHeight);
	glClear(GL_COLOR_BUFFER_BIT);
	installProgram(mapper);
	double phase = glutGet(GLUT_ELAPSED_TIME) * 5e-4;
	glUniform1f(glGetUniformLocation(mapper, "cos_phase"), cos(phase));
	glUniform1f(glGetUniformLocation(mapper, "sin_phase"), sin(phase));
	glUniform1i(glGetUniformLocation(mapper, "a"), 0);
	glUniform1i(glGetUniformLocation(mapper, "b"), 1);
	glUniform1i(glGetUniformLocation(mapper, "palette"), 2);
	glUniform1f(glGetUniformLocation(mapper, "fading"), fading);
	glActiveTextureARB(GL_TEXTURE1_ARB);
	glBindTexture(GL_TEXTURE_2D, transform[(currentTransform + 2) % 3]);
	glActiveTextureARB(GL_TEXTURE0_ARB);
	glBindTexture(GL_TEXTURE_2D, transform[(currentTransform + 1) % 3]);
	glActiveTextureARB(GL_TEXTURE2_ARB);
	glBindTexture(GL_TEXTURE_1D, palette);
	screen();
	glBindTexture(GL_TEXTURE_1D, 0);
	glActiveTextureARB(GL_TEXTURE0_ARB);

	// frame disaply
	installProgram(0);
	GLenum er = glGetError();
	if(er) printf("render: %s\n", ((char*)gluErrorString(er)));	
	glutSwapBuffers();
	glutPostRedisplay();
	
	// FPS counter
	static float delta = 0.1;
	delta = (chrono.time() + 9 * delta) * 0.1;
	if(F == 100) { F = 0;
	printf("%.2f FPS\n", 1.0 / delta);
	} else F++;
}

void reshapeFunc(int w, int h) {
	glViewport(0, 0, w, h);
	winWidth = w; winHeight = h; winRatio = w / (double)h;
}

int main(int argc, char** argv) {
	// open GL context setup //

	int err;
	glutInit(&argc, argv);

	glutInitWindowPosition(20, 20);
	glutInitWindowSize(768, 768);
	glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH | GLUT_STENCIL);
	glutCreateWindow("GPU fft");
   
	glutDisplayFunc(displayFunc);
	glutReshapeFunc(reshapeFunc);
	glutKeyboardFunc(key);
	//glutIdleFunc(renderScene);
   
	if(err = glewInit()) {
		printf("glew init: %s\n", glewGetErrorString(err));
		exit(1);
	}
	printf("OpenGL %s", glewGetString(GLEW_VERSION));
	if(GL_ARB_shading_language_100) {
		const char* v = (const char*) glGetString(GL_SHADING_LANGUAGE_VERSION_ARB);
		printf(" -- GLSL %s\n", v ? v : "1.0.51 (first version)");
	} else printf(" -- no GLSL\n");
	printf("frame buffer objects:\t%s\n", glGenFramebuffersEXT ? "yes" : "no");
	printf("point sprites:\t\t%s\n", GL_ARB_point_sprite ? "yes" : "no");
	printf("mipmap generation:\t%s\n", GLEW_SGIS_generate_mipmap ? "yes" : "no");
	GLint i;
	glGetIntegerv(GL_MAX_DRAW_BUFFERS, &i);
	printf("total bufferrs:\t%d\n", i);
	glGetIntegerv(GL_AUX_BUFFERS, &i);
	printf("aux bufferrs:\t%d\n", i);
	glGetIntegerv(GL_MAX_TEXTURE_UNITS_ARB, &i);
	printf("texture units:\t%d\n", i);
	
	#ifdef WIN32
	if(wglSwapInterval) wglSwapInterval(1);
	else if(wglSwapIntervalEXT) wglSwapIntervalEXT(1);
	#else
	//if(glXSwapInterval) glXSwapInterval(1);
	int (*SwapInterval)(int) = 0;
	SwapInterval = (int (*)(int)) glXGetProcAddress((GLubyte*)"glXSwapInterval");
	if (!SwapInterval)
	SwapInterval = (int (*)(int)) glXGetProcAddress((GLubyte*)"glXSwapIntervalEXT");
	if (!SwapInterval)
	SwapInterval = (int (*)(int)) glXGetProcAddress((GLubyte*)"glXSwapIntervalSGI");
	if(SwapInterval) SwapInterval(1);
	#endif
	
	
	// Shaders setup //
	mapper = createProgram(0,
		"uniform sampler2D a, b;"
		"uniform sampler1D palette;"
		"uniform float cos_phase, sin_phase, fading;"
		"void main() {"
		"vec2 complex = mix(texture2D(a, gl_TexCoord[0].st).rg, texture2D(b, gl_TexCoord[0].st).rg, fading);"
		"gl_FragColor = texture1D(palette, complex.r * cos_phase + complex.g * sin_phase);" 
		"}"
	);
	
	// fft setup //
	data = new float[4 * dataSize * dataSize];
	glGenTextures(3, transform);
	for(int i = 0; i < 3; i++) {
		glBindTexture(GL_TEXTURE_2D, transform[i]);
		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, dataSize, dataSize, 0, GL_RGBA, GL_FLOAT, 0);
		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
		glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
	}
	srand(Chrono::get_time());
	glGenTextures(1, &spectrum);
	glBindTexture(GL_TEXTURE_2D, spectrum);
	glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F_ARB, dataSize, dataSize, 0, GL_RGBA, GL_FLOAT, 0);
	textureSetup();
		
	fft = new FFT(dataSize, true);
	
	// init palette //
	const int paletteSize = 64;
	for(int i = 0; i < paletteSize * 4; i++) {
		double r = i & 1;
		double l = i / (double)paletteSize * (0.9 + 0.1 * cos(23 * (paletteSize * 4 - i) / (double)paletteSize));
		double s = 1 - l;
		double c[3];
		double theta = i / (double)paletteSize * 2;
		double phi = i / (double)paletteSize * 5 + 1;
		c[0] = cos(theta) * cos(phi)  + 0.1 * rand() / (double)RAND_MAX;
		c[2] = sin(theta) * cos(phi)  + 0.1 * rand() / (double)RAND_MAX;
		c[2] = sin(phi) + 0.1 * rand() / (double)RAND_MAX;
		double n = 0;
		for(int j = 0; j < 3; j++) n += c[j] * c[j];
		for(int j = 0; j < 3; j++) data[i * 4 + j] = l * (s * c[j] / n + 1 - s);
		data[i * 4 + 3] = 0;
	}
	glGenTextures(1, &palette);
	glBindTexture(GL_TEXTURE_1D, palette);
	glTexImage1D(GL_TEXTURE_1D, 0, GL_RGBA, paletteSize, 0, GL_RGBA, GL_FLOAT, data);
	glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
	glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
	glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
	glTexParameterf(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
	
	
	// Main loop control //
	GLenum er = glGetError();
	if(er) printf("init: %s\n", (char*)gluErrorString(er));
	puts("--- starting ---");
	glutMainLoop();
   
}
